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1. ABSTRACT 

The feasibility of a viscoplastic model incorporating two back stresses and a 
drag strength is investigated for performing nonlinear finite element analyses of 
structural engineering problems. The model has been recently put forth by Freed 
and Walker. To demonstrate suitability for nonlinear structural analyses, the 
model is implemented into a finite element program and analyses for several 
uniaxial and multiaxial problems are performed. Good agreement is shown between 
the results obtained using the finite element implementation and those obtained 
experimentally. The advantages of using advanced viscoplastic models for 
performing nonlinear finite element analyses of structural components are indicated. 

2. INTRODUCTION 

Classical creep-plasticity constitutive models treat creep and plastic strains 
as independent noninteracting entities. These models are, therefore, incapable of 
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accounting for the observed interactions between creep and plastic strains at high 
temperatures. Viscoplastic models, however, consider all the inelastic strain 
(including plasticity, creep, relaxation, etc.) as a single, unified, time-dependent 
quantity, and thus, automatically include interactions that occur among them. 
Viscoplastic models, therefore, provide more realistic descriptions of time- 
dependent inelastic behavior of materials at high temperatures. 

Viscoplastic models become more realistic when as much material physics 
as possible is included in the model. This, however, results in complex 
mathematical frameworks for viscoplastic models. The constitutive differential 
equations of viscoplastic models that govern the flow and evolution laws are 
generally highly nonlinear and mathematically “stiff.” The closed-form solutions 
for structural engineering problems are virtually intractable when viscoplastic 
models are used to define the stress-strain relationship. To assess the advantages 
offered by more realistic viscoplastic models one must, therefore, employ 
numerical solution methodology involving, for example, the finite element method 
or the boundary element method. 

Toward this aim, this paper presents a finite element solution methodology 
developed at the National Aeronautics and Space Administration (NASA) Lewis 
Research Center for use with the general purpose finite element program MARC 
[MARC Analysis Research Corporation, 1992]. The methodology, designed for 
use with viscoplastic models, is demonstrated in this paper for a viscoplastic 
model recently developed by Freed and Walker [1993] that incorporates two back 
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stresses. However, because the methodology is general in nature, it can easily be 
adapted for use with other viscoplastic models. 

For completeness, the paper includes brief descriptions of the viscoplastic 
model by Freed and Walker [1993] and the finite element solution methodology. 
The methodology is illustrated by applying it to several uniaxial and multiaxial 
problems. 

3. VISCOPLASTIC MODEL 

The viscoplastic model used herein was recently put forth by Freed and 
Walker [1993]. The model contains one scalar internal state variable D, called the 
drag strength, and a tensorial internal variable B- , called the back stress. The 
back stress B- is assumed to be composed of two back stresses that are denoted 
by B®j and B?j. (A small displacement and a small strain formulation is 
employed in the model.) 

The stress a - is taken to be related to elastic strain by Hooke’s law 


<; = L±^ij - *ij (iJ = 1,2,3) 


where E is Young’s modulus, v is Poisson’s ratio, and is the Kronecker delta. 
The symbol e denotes the strain and superscript e denotes the elastic component 
of strain. Following Einstein’s summation convention, the repeated subscripts in 
equation (l) and elsewhere imply summation over the range of the symbol. 

The total strain rate £•• is written as the sum of elastic £?•, inelastic 

V J 

(including plasticity, creep, relaxation, etc.) eY., and thermal £- strain rate 
components. In symbol form, 


3 



*U = *'j + + *‘ U (i,j = 1,2,3) 


In equation (2) the superscripts v and t denote the viscoplastic and thermal 
components of strain rate, respectively. A dot over a symbol indicates its 
derivative with respect to time t. 

The deviatoric stress S— and the deviatoric total strain rate E- have the 
following expressions 


S ij a ij 3 ff kk ^ij 


1 

E ij = j ^kk *ij 

The back stress B-j is the sum of two back stress components B- and B-, 


that is 


B:: = B S . + B*. 

1J »J 


The effective stress S - is defined as 


(6) 2.. = S-- - B 

The temperature dependence of the model is mainly contained in the 
thermal diffusivity function 0, which is defined as 


exp(— Q/kT) 


T t < T < T m 


T t 

In — +1 


0 < T <T t 
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In equation (7), Q is the activation energy, k is the Boltzmann constant, T 
is absolute temperature, T t is transition temperature, and T m is the melting 
temperature of the material. 

The function Z, the Zener- Hollomon parameter, is defined by 


( 8 ) 


Z = A sinh 11 



where 

h = s ij V 2 

Here A and n are inelastic material constants. 

The evolution equations for the internal state variables are 


( 9 ) 


B . = 2H 

ij c 


- 2L 


B“. 

— Io 


a 


where a = s,{. The strain invariant I 2 is defined by 


( 10 ) 


I 2 = 1 


ez 


[ P ijV P 


M 


In equation (10) 


( 11 ) 


P : : = 


E 


if J 2 < 

otherwise 


S ij 


U 2(1 +v) 


and E” is defined by equation (4). Also 


( 12 ) 


P M - 


E 


2(1 + i/) 




K - - B“ Sij 


,j ^n 


The functions K, H ft , and in the above equations are given by 


o 


(13) 


K = (C - D 0 ) (C + D 0 )/4£C 
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(14) 
and 

(15) 


H s = 


E 


2(1 + i/) 


H« = 


E 


2(1 + u) H 0 


L s = f 


(C - D) (D - Dq) 


h = (1 " 3 


(C - D) (D - D 0 ) 


in which C, 6, H 0 , and f are material constants and D 0 denotes the minimum 
value that the drag strength D can take for a given material. 

The drag strength D evolves according to the following equation 


(16) 

where 


D= h(I 2 - ©r) D 0 < D < D 


max 


= (C + D 0 )/2 


and 


r = A sinh [(D - D 0 )/£C] n 
The function h is defined as 


(17) 


where 


(18) 


h = h 


D 


(D - D 0 )/6C 


sinh{(D - D 0 )/^C} 


h D - 


ho - — T 


T t < T <T m 


0 < T < T t 


In equation (18), h Q and hj are material constants. 

In addition to elastic material constants, the shear modulus /x, Poisson’s 
ratio v , and the coefficient of thermal expansion a, the above viscoplastic model 
has 13 inelastic material constants: A, C, £, D 0 , f, H 0 , h Q , h 1? m, n, Q, T , and 
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T f . The values of these constants taken from Freed and Walker [1993] are listed 
in Table I. 

4. DESCRIPTION OF MARC PROGRAM 
A viscoplastic model is most advantageous if it can be used on a practical 
scale for the solution of structural engineering problems faced by, for example, 
the aerospace and nuclear industries. To demonstrate the feasibility of the 
present viscoplastic model for performing complex nonlinear structural analyses, 
the model was implemented into the general purpose finite element program 
MARC [1992]. Any other finite element program capable of performing nonlinear 
viscoplastic structural analyses can also be used. The MARC program employs 
sophisticated integration algorithms and advanced finite element formulations. It 
is specially tailored to fit the requirements of a nonlinear structural analysis. 

The MARC program provides the user with a convenient way of 
implementing and integrating the constitutive differential equations of F reed and 
Walker’s model [F & W, 1993]. This is accomplished through the user subroutine 
HYPELA [see MARC Analysis Research Corporation, 1992] provided in MARC. 
The subroutine HYPELA, developed for Freed and Walker’s model, as well as 
the self-adaptive time integration strategy used in it, follow essentially the same 
structure as given in Ary a and Kaufman [1989] and Cassenti [1983]. 

In the MARC program, all the material nonlinearities are incorporated into 
an initial load vector that is treated as a pseudo-body force in finite element 
equilibrium equations. The constitutive equations of viscoplastic models are 
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mathematically “stiff.” A subincrement technique is employed to form the 
incremental constitutive equations corresponding to a given finite load increment. 

In the subincrement technique, a given finite load increment is split into 
several equal sub increments. The constitutive equations of the viscoplastic models 
are then integrated over these small subincrements to obtain an accurate 
representation of the incremental constitutive equations over the finite load 
increment. Structural engineering problems encountered in aerospace and nuclear 
industries usually involve many degrees of freedom. Since the use of an explicit 
integration method does not require the assembly and inversion of Jacobian 
matrices, it is more suitable for problems involving a large number of degrees of 
freedom. The explicit forward Euler method with a self-adaptive time step was 
used to integrate the constitutive equations over the subincrements. The 
integration strategy is found to work efficiently and accurately, even for large 
finite element load increments, provided the subincrements are suitably small to 
ensure the stability of the forward Euler method (see, for example, Arya [1992] 
and Arya and Arnold [1992]). However, it is difficult for the user to select 
efficient subincremental steps, and there is a considerable incentive (saved CPU 
time) to use as few subincrements as possible without affecting the stability of 
the constitutive equations associated with the viscoplastic model. For an 
experienced user, however, it is easier to optimize the number of subincrements. 

A brief introduction to the MARC program is provided here to familiarize 
the user with its operation. (For complete details of the MARC program, see 
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MARC Analysis Research Corporation, 1992.) This introduction, taken from 
Walker [1981], is as follows: 

“The principle of virtual work may be used to generate the 
MARC non-linear equilibrium equations governing the incremental 
response of the structure to an increment of load. In evaluating the 
non-linear structural response of a component, the program assumes 
that the load history is divided into incrementally applied loads. 

Each load step is sequentially analyzed as a linear matrix problem 
using an appropriate stiffness matrix and load vector. Although 
each load step uses linear matrix methods to solve the incremental 
equations, the incremental equations themselves are non-lineax since 
the load vector will depend on the displacement increment obtained 
in the solution of incremental equilibrium equations. 

The principle of virtual work may be written, for applied 
external point loads Pj, or displacement Uj, in the form: 

E f ££ T <T: dv = <Su T P: (^-l) 

JV 1 1 1 

where the integral extends over the volume, V, of each finite 
element and the summation sign extends to all the elements in 
structure. 

In Eq. (5.1) the virtual displacement vector £uj is related to 
the virtual strain vector through the relationship: 
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6 6: = B:{6 n j so that tie? = ^.bT. (5*2) 

1 j j u j i uj ij 

where B- is the strain displacement matrix and the superscript T 

denotes transposition. Since <5u- is an arbitrary virtual displacement 

vector, Eqs. (5.1) and (5.2) may be written in the form: 

I f B T .<7: dV = P: (5.3) 

Jv IJ J 1 

This relation expresses the equilibrium of structure when the 
applied load vector is P- and the stress vector is a- v If an 
incremental load APj is applied to the structure and the stress 
vector changes to cr- + A a-, the relation expressing the equilibrium 
of the structure at the end of the incremental load application may 
be written as: 

1 Jv B U (<, i + A,7 i> dv = p i + AP i (5 ' 4) 

Hence, the relation expressing the equilibrium of the structure 
during the application of the incremental load vector APj is 
obtained from Eqs. (5.3) and (5.4) by subtraction in the form: 

I f B T . A<7; dV = AP: (5.5) 

J V iJ J 1 

The MARC code allows the user to implement very general 
constitutive relationships into the program by means of the user 
subroutine HYPELA. Within this subroutine, the user must specify 
the values of elasticity matrix D- and the inelastic stress vector A{; 
in the incremental vector constitutive relationship: 
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A<7j = DjjfAej - S^a A e) - ACi ( 5 * 6 ) 

The inelastic stress increment vector ACj is computed in HYPELA 

using the constitutive relationships ... [of the viscoplastic model] . 

In Eq. (5.6), a denotes the coefficient of thermal expansion 

and is the vector Kronecker delta symbol, 

C f 1 if 0 < j < 3 (5.7) 

J \ 0 if 3 < j < 6 

For the . . . [viscoplastic models], the incremental inelastic 
stress vector A£j depends in a highly nonlinear manner on the 
incremental strain vector Aej. Since Asj = B”A u j, the incremental 
inelastic stress vector AC- depends in a highly non-linear manner 
on the nodal displacement vector A u j, so that ACj = AC(A u j). 

Substitution of Eq. (5.6) into (5.5) produces the incremental 
equilibrium equations for MARC program in the form: 

EK ij Au j = AP i + AR i +S JyByACjdV + S J v B^<jaA»dV 

(5.8) 

where Ky is the elemental elastic stiffness matrix defined by the 
relation: 

K U = Iv B fk D k*B,j dV (5.9) 

The vector AR- is the residual load correction vector or out-of- 
equilibrium force vector from the preceding load increment: 

ARj = Pj - E J y B^j dV (5.10) 
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which is added to the current increment in order to restore the 
structure to equilibrium. The nonlinearity in the incremental 
equilibrium relationship, defined in Eq. (5.8), arises because the 
inelastic stress increment vector A£- depends non-linearly on the 
displacement increment vector Auj. Values of D-j and ACj 
appropriate to the current incremental load step are returned to the 
main program by subroutine HYPELA and the incremental 
equilibrium relations in Eq. (5.8) are solved by successive iterations. 

The solution of the incremental equilibrium equations in (5.8) 
is accomplished within the MARC code by the following algorithm. 
At the start of the increment the user subroutine HYPELA is 
entered to determine the elasticity matrix D- and the incremental 
inelastic stress vector A£j. On entry to the subroutine, the input 
consists of the strain increment vector Asj, the temperature 
increment A©, the time increment At over which the incremental 
external load vector APj is applied to the structure, and the values 
of stress, strain, temperature and viscoplastic state variables at the 
beginning of the increment. Since the incremental strain vector 
Asj = BjjA u j, can only be accurately determined after the solution 
to the incremental equilibrium relationship in Eq. (5.8) has yielded 
the correct incremental solution vector A u j, the strain increment 
vector Aej initially used to generate the inelastic stress increment 
vector A^j must be estimated. The initial estimate for A£j is 
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assumed to be the value obtained for Ae- in the preceding incre- 
ment. On exit from subroutine HYPELA the elasticity matrix D- 
and the estimated inelastic stress increment vector ACj are passed 
into the main program. After the values of D- and A £ j are 
obtained for each integration point in the structure, the incremental 
equilibrium relationship in Eq. (5.8) is assembled and solved for the 
incremental node displacement vector A u j. The incremental strain 
vector, = By A u j, is then computed and compared with the 
initial guess for Ae- used to generate the inelastic incremental stress 
vector ACj. If this incremental strain vector is equal, within a user 
specified tolerance, to the incremental strain vector used to 
compute ACj in the assembly phase, the solution is assumed to 
have converged. Otherwise the updated strain increment vector, 
obtained from the solution of the equilibrium relations in Eq. (5.8), 
is passed into subroutine HYPELA, a new vector, ACj, is computed 
and the equilibrium equations resolved to yield an improved value 
of A u j and Ae^. The process is repeated until the value of vector 
AfiN on the assembly phase is equal, within a user specified 
tolerance, to the value of vector Ae^ on the solution phase. After 
convergence is achieved, the temperature, stress vector, strain 
vector and viscoplastic state variables are updated by adding the 
incremental values generated during the current increment to the 
values of these variables at the beginning of the increment. The 
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program then passes on to the next load increment where the 
process is repeated.” 

5. APPLICATION TO PROBLEMS 

5.1 Uniaxial Problems 

Some uniaxial problems were solved first to demonstrate and validate the 
implementation of Freed and Walker’s model into the MARC program. These 
included calculations of tensile stress-strain response at different temperatures 
and cyclic response under isothermal and nonisothermal loadings. The 
constitutive equations (including the flow and evolution laws) of the viscoplastic 
model were implemented in MARC in their generalized three-dimensional form 
through the user subroutine HYP EL A. This generalized finite element 
implementation was used to perform the uniaxial computations that are 
presented in this paper. The values of constants for copper listed in Table I and 
taken from Freed and Walker [1993], were utilized in these computations. 

Figures 1 through 4 show the results for uniaxial problems. 

5.2 Multiaxial Cowl Lip Problem 

To demonstrate the feasibility of performing complex structural analyses 
using a viscoplastic model, the validated finite element implementation of Freed 
and Walker’s viscoplastic model was applied to a sample structural component 
used in the aerospace industry. The component, called a cowl lip, is part of the 
leading edge of a hypersonic aircraft engine inlet. To achieve high inlet 
aerodynamic performance for proposed hypersonic flight between Mach 3 and 25, 
not only must the high heat flux and high heating rates be tolerated, but 
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distortions caused by thermal warping of the structure must also be minimized. 
Consequently, the need arises for the development of actively cooled leading 
edges fabricated from specialized materials as well as innovative cooling concepts 
that enable a structure to withstand the severe thermal loading conditions. The 
details of different cooling concepts proposed under a NASA Lewis sponsored 
program, called COLT (Cowl Lip Technology Program), were presented by Melis 
and Gladden [1988]. In this paper, a cooling concept, called the parallel flow 
concept, is investigated. In this concept, the coolant channels are laid parallel to 
the leading edge of the cowl lip (fig. 5), and the coolant is flowed through them 
to contain the temperature of the component. The viability of the parallel flow 
concept is investigated by performing the structural analysis with the finite 
element implementation of Freed and Walker’s viscoplastic model. 

The cowl lip geometry is shown in figure 5. The dimensions of the cowl lip 
are 15.2 by 3.8 by 0.64 cm. However, for the finite element analysis only the 5-cm 
central portion of the cowl lip was modeled to avoid end effects that would be 
difficult to quantify. The finite element mesh for this central portion was 
constructed by Melis and Gladden [1988]. The model was made up of 3294 solid, 
8-noded hexagonal elements and has 4760 nodes. A large number of finite 
elements was required to deal with the severity of thermal loadings and 
gradients. Note from figure 6 that the temperature at a critical location of the 
component rises from 21 to 758 °C in only 0.75 sec, indicating the severity of 
loading. However, because of the geometrical symmetry of the component, one 
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need consider only half of the component for the finite element analysis, which 
results in a considerable savings in CPU time. 

The steady-state temperature distribution in the component was obtained 
from Melis and Gladden [1988], which includes detailed results of a steady-state 
heat transfer analysis that was performed to obtain the steady-state temperature 
distribution in the cowl lip. Figure 6 shows this temperature distribution in the 
cowl lip. Figure 7 depicts the thermal loading cycle used in finite element 
analysis. The transient temperature distribution was obtained by using a linear 
interpolation technique. The highest temperatures in the component occur along 
the leading edge of the cowl lip. The temperature values thus obtained for the 
complete thermal loading cycle were used to perform the cyclic finite element 
analysis of the cowl lip. 

6. RESULTS AND DISCUSSION 

6.1 Uniaxial Problems 

Figure 1 shows isothermal tensile stress-strain curves at different 
temperatures obtained by using the finite element implementation of Freed and 
Walker’s viscoplastic model. The stress rate used in these computations is 
0.0061 MPa s' 1 . Figure 1 also shows a comparison of curves obtained by Freed 
and Walker with those obtained experimentally. Good agreement exists among 
the curves obtained by Freed and Walker, by experiments, and by using the 
present implementation of the viscoplastic model. 

The thermal and mechanical loadings used to generate the isothermal and 
nonisothermal hysteresis loops are shown in figure 2. These hysteresis loops for 
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isothermal and nonisothermal loadings are shown in figures 3 and 4, respectively. 
The strain rates used in generating these isothermal and nonisothermal hysteresis 
loops are 0.001 s' 1 and 1.5 xlO’ 5 s’ 1 , respectively. These figures also exhibit the 
hysteresis loops obtained by experiments. Again good agreement exists between 
the curves obtained experimentally and by using the finite element 
implementation. The difference in the values of stress obtained by using the finite 
element implementation and those observed in experiments (fig. 4), is due to the 
viscoplastic model. An error in the finite element implementation is not implied 
because the present results are in excellent agreement with the results of Freed 
and Walker [1993]. 

The good agreement between the calculated values from the finite element 
implementation and those obtained from experiments for the tensile curves and 
isothermal and nonisothermal hysteresis loops validates the finite element 
implementation of Freed and Walker’s model. Such validation is a worthwhile 
investment of effort and computer time before applying the model to industrial 
structural engineering problems. These problems, in general, involve complex 
thermal/mechanical loadings in addition to their complex geometries and are, 
therefore, computationally intensive and time consuming. 

6.2 Cowl Lip Problem 

The most significant results of the cowl lip problem were obtained by 
plotting the stresses and strains along the leading edge, that is, the z-direction in 
figures 8 through 11. The stress distribution in the z-direction at 0.75 s is 
exhibited in figure 8, which shows the maximum (compressive) stress occurring 
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at the leading edge of the cowl lip. Figure 9 displays the stress distribution in the 
z-direction of the cowl lip at 2.25 s. A comparison of these figures reveals 
significant redistribution of stress in the component as a result of inelastic 
deformation. The stress along the leading edge is now seen to be tensile. This 
shows that viscoplastic models are capable of picking the redistribution of stress 
caused by inelastic deformation even for a short duration of 1.5 s. Note that a 
nonunified elastic-plastic-creep analysis of this component performed by Arya 
et al. [1991] was unable to capture this redistribution of stress. 

The strain distributions in z-direction in the cowl lip at 0.75 s and 2.25 s 
are shown in figures 10 and 11, respectively. These figures show the values of 
total strain at different locations of the cowl lip. The compressive inelastic strain 
along the leading edge accumulates (increases) with time, which makes the total 
(tensile) strain along the leading edge decrease. This can be observed by 
comparing figures 10 and 11. For example, the total strain along the edge of the 
cowl lip reduces from a value of 0.00491 at 0.75 s to a value of 0.00483 at 2.25 s. 

The deformed shape of the segment at 0.75 s is plotted in figure 12. For 
easier examination, the displacements in this figure are magnified by a factor of 
1000. This figure shows that the maximum deformation of the segment occurs 
near point A. The component thickens in this region, which will lead to 
distortion of the coolant channel configuration. This thickening of the leading 
edge will also require more cooling down time, and this will lead to an early 
failure of the component. This indicates that, for the cowl lip problem, the 
parallel flow concept is not a particularly advantageous concept. 
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7. CONCLUSIONS 


A n um erical solution methodology based on the finite element method is 
described for a viscoplastic model recently developed by Freed and Walker. It is 
illustrated for the finite element program MARC. However, the generality of 
methodology makes it easily adaptable to any other finite element program. The 
viability of methodology is demonstrated by applying it to several uniaxial 
problems and a multiaxial cowl lip problem. The results for uniaxial tensile and 
cyclic loadings obtained by using the finite element implementation show good 
agreement with the experimental results and with the results reported in the 
literature. The advantages of using an advanced viscoplastic model for nonlinear 
structural analyses are established by applying it to a multiaxial cowl lip problem 
and investigating the proposed parallel flow design concept. The results from the 
viscoplastic finite element analysis indicate that the parallel flow concept for the 
design of cowl lip is not a particularly advantageous concept, as it will lead to an 
early failure of the component. 
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TABLE I.— MATERIAL CON- 
STANTS FOR COPPER 
[F&W, 1993] 


[n = Hq + /ijT; T is in degrees 
Kelvin; D Q = C/100; S = 0.035; 
T, = 0.5 T.] 

£ mJ 


Constant 

Unit 

Value 

a 

K' 1 

18xl0’ 6 

M 0 

MPa 

43 000 


MPa/K 

-17 

V 

— 

0.36 

A 

s" 1 

2xl0 7 

C 

MPa 

13 

f 

— 

0.75 

H 0 

— 

20 

^0 

MPa 

50 


MPa 

15 

m 

— 

0.5 

n 

— 

4.5 

Q 

J/mol 

200 000 

T 

m 

K 

1356 
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Strain 

Figure 1.— Tensile stress-strain curves at differ- 
ent temperatures. (Stress rate is 0.0061 
MPa/s.) 
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Figure 3.— Isothermal hysteresis loops at 21 °C. 
(Strain rate is 0.02/s; strain range is 4 percent.) 
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Figure 2.— Thermal and mechanical loadings used to generate 
hystersis loops by using the finite element implementation. 
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Figure 4. — Nonisothermal (in phase) hysteresis 
loops. (Temperature range is 200 to 500 °C; 
strain range is 0.0075; and strain rate is 
1.5 x 10-5/s.) 
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Figure 8. — Stress in the cowl lip at 0.75 s; z-component. 
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Figure 5. — Cowl Hp finite element model (3294 elements; 4760 
nodes). 
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Figure 9. — Stress in the cowl lip at 2.25 s; z-component. 
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Figure 6. — Steady-state temperature distribution in the cowl lip. * 


Figure 1 0. — Strain in the cowl lip at 0.75 s; z-component. 
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Figure 7. — Thermal loading cycle. 
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Figure 11. — Strain in the cowl lip at 2.25 s; z-component. 



Figure 1 2.— Deformed shape of the segment at 0.75 s. (For easier 
examination, this deformation is magnified by a factor of 1000.) 
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